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Abstract 


This  study  describes  different  analytical  models  that  have  previously  been  developed  for 
predicting  the  radial  growth  and  collapse  of  underwater  explosion  (UNDEX)  gas  bubbles 
in  a  free-field  environment.  The  report  describes  the  implementation  of  nine  analytical 
gas  bubble  models,  in  the  form  of  nonlinear  differential  equations,  and  a  fourth-order 
Runge-Kutta  solution  method.  Gas  bubble  radius  time  histories  calculated  with  these 
models  are  compared  to  empirical  models  derived  from  published  experimental  data.  The 
analytical  models  allow  for  different  assumptions  such  as  fluid  compressibility,  bubble 
migration  coupled  to  dilatation,  and  an  empirical  correction  for  energy  loss.  It  was  found 
that  none  of  the  analytical  models  fully  account  for  the  reduction  in  the  gas  bubble  radius 
throughout  the  growth  and  collapse  cycles.  Including  compressibility  in  the  fluid  and  the 
gas  bubble  provides  the  best  predictions  when  compared  to  experimental  fits.  The 
incompressible  fluid  model  requires  an  empirical  energy  loss  function,  as  there  is  no 
energy  loss  inherent  within  the  model.  Models  considering  just  the  compressibility  of  the 
surrounding  fluid  do  not  account  for  the  full  energy  loss  seen  in  the  experimental  fits,  and 
produced  similar  results.  Inclusion  of  migration  effects  had  no  influence  on  the  bubble 
radius  or  period  because  of  the  large  detonation  depth. 


Resume 


La  presente  etude  decrit  differents  modeles  analytiques  elabores  anterieurement  pour 
prevoir  la  croissance  et  T  implosion  radiales  en  champ  fibre  des  bulles  gazeuses  creees  par 
une  explosion  sous-marine  (UNDEX).  Le  rapport  decrit  la  mise  en  oeuvre  de  neuf 
modeles  analytiques  de  bulles  gazeuses  utilisant  des  equations  differentielles  non 
lineaires  et  une  methode  de  solution  de  Runge-Kutta  du  quatrieme  ordre.  On  compare 
T  evolution  temporelle  du  rayon  des  bulles  gazeuses  calculee  avec  ces  divers  modeles 
avec  celle  obtenue  avec  des  modeles  empiriques  derives  de  donnees  experimentales 
publiees.  Les  modeles  analytiques  pennettent  differentes  hypotheses  concernant,  entre 
autres,  la  compressibilite  du  fluide,  la  migration  des  bulles  reliee  a  la  dilatation,  et  la 
correction  empirique  de  la  perte  d’energie.  On  constate  qu’aucun  des  modeles  analytiques 
ne  peut  expliquer  completement  la  reduction  du  rayon  des  bulles  gazeuses  sur  toute 
Tetendue  des  cycles  de  croissance  et  d’implosion.  L’inclusion  des  compressibilites  du 
fluide  et  des  bulles  gazeuses  donne  les  meilleures  previsions  comparativement  aux 
ajustements  avec  les  donnees  experimentales.  Le  rnodele  a  fluide  incompressible  exige 
d’utiliser  une  fonction  empirique  de  perte  d’energie  etant  donne  qu’il  ne  tient  pas  compte 
lui-meme  de  la  perte  d’energie.  Les  modeles  qui  ne  tiennent  compte  que  de  la 
compressibilite  du  fluide  environnant  ne  peuvent  expliquer  completement  la  perte 
d’energie  observee  dans  les  ajustements  aux  donnees  experimentales  et  donnent  des 
resultats  similaires.  L’inclusion  des  effets  dus  a  la  migration  n’a  aucun  effet  sur  le  rayon 
ou  la  duree  des  bulles  a  cause  de  la  grande  profondeur  a  laquelle  les  explosions  ont  lieu. 
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Executive  summary 


Analytical  Solutions  for  Predicting  Underwater  Explosion  Gas 
Bubble  Behaviour 

Mark  Riley;  DRDC  Atlantic  TM  2010-237;  Defence  R&D  Canada  -  Atlantic; 

November  2010. 

Introduction:  Ship  vulnerability  to  underwater  explosion  loading  is  being  investigated 
as  part  of  the  Maritime  Force  Protection  (MFP)  TDP.  The  ultimate  goal  of  this  project  is 
to  deliver  a  validated  capability  to  predict  these  types  of  loading  which  can  then  be 
applied  to  full  scale  ships.  This  study  developed  a  computer  code  to  implement  analytical 
models  that  will  provide  a  quick  and  easy  way  to  determine  the  radial  and  migration 
behaviour  of  an  underwater  explosion  gas  bubble.  The  analytical  results  are  compared  to 
experimental  fits. 

Results:  This  report  compares  various  analytical  models  for  predicting  the  growth  and 
collapse  of  an  underwater  explosion  gas  bubble.  The  different  models  use  various 
underlying  assumptions  such  as  compressibility  of  the  surrounding  fluid  and  gas  bubble, 
energy  losses,  and  the  coupling  of  radial  and  migration  motion.  This  study  shows  that  the 
reduction  in  gas  bubble  radius  through  the  growth  and  collapse  cycles  is  underestimated 
by  all  of  the  analytical  models  when  compared  to  experiments.  Of  the  models  trialed,  the 
best  results  were  obtained  when  compressibility  of  the  gas  and  surrounding  fluid  was 
included. 

Significance:  The  analytical  models  trialed  in  this  study  are  simple  alternatives  to  two 
dimensional  and  three-dimensional  numerical  simulations  involving  run  times  up  to  14 
days  on  30-40  processors.  Using  analytical  models,  the  run  time  is  in  the  order  of 
seconds  for  the  same  simulation  time.  This  study  is  an  initial  step  in  the  development  of 
a  code  that  will  predict  the  loading  from  an  underwater  explosion  gas  bubble  in  close 
proximity  to  a  target.  The  results  of  this  study  show  that  there  are  various  underlying 
assumptions  that  can  be  made  which  significantly  affect  the  gas  bubble  response. 

Future  plans:  The  models  implemented  in  this  study  could  be  extended  to  predict  far- 
field  bubble  loads  on  ships  and  subs  for  analysis  such  as  whipping.  Further  development 
of  the  code  is  planned  to  include  the  influence  of  a  close  proximity  target  on  the 
behaviour  of  the  gas  bubble. 
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Analytical  Solutions  for  Predicting  Underwater  Explosion  Gas 
Bubble  Behaviour 

Mark  Riley;  DRDC  Atlantic  TM  2010-237;  R  &  D  pour  la  defense  Canada  - 

Atlantique;  Novembre  2010. 

Introduction  :  La  vulnerabilite  des  navires  a  la  charge  exercee  par  les  explosions  sous- 
marines  est  actuellement  etudiee  dans  le  cadre  du  Programme  de  demonstration  de 
technologies  (PDT)  de  protection  de  la  force  navale  (PFN).  Le  but  ultime  du  projet  est  de 
mettre  au  point  un  outil  valide  permettant  de  calculer  les  charges  de  ce  type,  puis  de 
l’appliquer  a  des  navires  de  grandeur  nature.  La  presente  etude  a  elabore  un  code  machine 
pour  mettre  en  oeuvre  des  modeles  analytiques  qui  permettront  de  determiner  rapidement 
et  facilement  les  comportements  radial  et  migratoire  d’une  bulle  gazeuse  creee  par  une 
explosion  sous-marine.  On  compare  les  resultats  analytiques  avec  les  ajustements  aux 
donnees  experimentales. 

Resultats  :  Le  present  rapport  compare  divers  modeles  analytiques  servant  a  prevoir  la 
croissance  et  l’implosion  d’une  bulle  gazeuse  creee  par  une  explosion  sous-marine.  Ces 
modeles  font  diverses  hypotheses  concernant,  entre  autres,  les  compressibilites  du  fluide 
environnant  et  de  la  bulle  gazeuse,  la  perte  d’energie,  et  le  couplage  des  mouvements 
radial  et  migratoire.  L’etude  montre  que,  comparativement  aux  resultats  d’experience,  la 
reduction  du  rayon  de  la  bulle  gazeuse  dans  les  cycles  de  croissance  et  d’implosion  est 
sous-estimee  par  tous  les  modeles  analytiques.  Parmi  les  modeles  mis  a  l’essai,  on  a 
obtenu  les  meilleurs  resultats  quand  on  a  tenu  compte  des  compressibilites  du  gaz  et  du 
fluide  environnant. 

Importance  :  Les  modeles  analytiques  mis  a  l’essai  dans  la  presente  etude  sont  des 
modeles  simples  visant  a  remplacer  des  simulations  numeriques  bi-  et  tridimensionnelles 
d’une  duree  d’execution  pouvant  atteindre  14  jours  sur  30  a  40  processeurs.  Avec  ces 
modeles  analytiques,  la  duree  d’execution  est  de  l’ordre  de  quelques  secondes.  L’etude 
est  un  premier  pas  dans  l’elaboration  d’un  code  qui  predira  la  charge  exercee  par  une 
bulle  de  gaz  creee  par  une  explosion  sous-marine  au  voisinage  immediat  d’une  cible.  Les 
resultats  de  l’etude  montrent  que  l’on  peut  s’appuyer  sur  diverses  hypotheses  et  que 
celles-ci  ont  un  effet  important  sur  le  comportement  de  la  bulle  gazeuse. 

Perspectives  :  Les  modeles  elabores  dans  la  presente  etude  pourraient  etre  appliques  a  la 
prevision  des  charges  des  bulles  en  champ  lointain  sur  les  navires  et  les  sous-marins  en 
vue  d’analyser  des  effets  comme  le  fouettement.  On  prevoit  ameliorer  le  code  pour  y 
inclure  1’ effet  d’une  cible  tres  proche  sur  le  comportement  de  la  bulle  gazeuse. 
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1  Introduction 


Naval  vessels  can  be  subjected  to  severe  blast  loading  from  both  air  and  underwater 
explosions.  For  underwater  explosions  in  close  proximity  to  a  structure,  the  loading 
mechanism  is  entirely  different  than  for  air  blast.  With  underwater  explosions,  a  shock  is 
fonned  that  propagates  through  the  water  and  loads  the  structure  with  a  significant 
pressure  pulse.  An  underwater  explosion  will  also  generate  a  pulsating  gas  bubble  whose 
size  depends  on  the  chemical  composition  of  the  explosive,  the  mass  of  the  charge,  and 
the  depth  of  the  charge  in  the  water.  When  the  charge  is  detonated  within  2-3  times  the 
maximum  free-field  bubble  radius,  R,  of  the  structure,  the  bubble  collapse  and  the  water 
jet  that  is  created  by  the  collapse  are  very  likely  to  be  directed  at  the  structure.  At  even 
closer  ranges  the  bubble  collapses  directly  onto  the  structure,  resulting  in  loading  that  can 
be  more  severe  than  the  shock  loading.  The  bubble  can  then  reform  and  pulsate  against 
the  structure,  resulting  in  more  loading  cycles.  During  the  bubble  collapse  the  boundary 
of  the  bubble  farthest  from  the  structure  typically  collapses  inward,  allowing  a  water  jet 
directed  toward  the  structure  to  fonn  and  causing  the  bubble  to  change  from  a  spherical  to 
toroidal  shape  as  shown  in  Figure  1. 


Figure  1:  Bubble  collapse  and  water  jet  formation  against  a  structure 

As  part  of  the  Maritime  Force  Protection  (MFP)  TDP,  ship  vulnerability  to  underwater 
explosion  loading  is  being  investigated.  The  ultimate  goal  of  this  project  is  to  deliver  a 
validated  capability  to  predict  these  types  of  loading  which  can  then  be  applied  to  full 
scale  ships.  Although  numerical  analysis  tools  are  available  for  modelling  the  loading 
and  structural  response  for  close-proximity  UNDEX,  these  methods  are  laborious  and 
very  time  consuming  and  still  exhibit  quantitative  gaps  when  compared  to  experimental 
results.  This  study  is  an  initial  step  in  developing  a  simplified  approach  to  determining 
the  loading  on  a  structure  due  to  close  proximity  UNDEX  gas  bubble  collapse. 
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This  report  examines  various  analytical  models  that  are  available  for  predicting  the  gas 
bubble  growth  and  collapse.  The  similitude  equations  for  explosive  materials  have  been 
shown  to  provide  good  results  for  various  shock  parameters,  such  as  the  peak  pressure, 
decay  time,  and  energy  as  well  as  the  gas  bubble  period  and  maximum  free-field  gas 
bubble  radius.  Major  limitations  of  the  similitude  equations  and  current  analytical 
approaches  are  they  do  not  include  effects  of  loading  from  the  gas  bubble  collapse,  they 
do  not  account  for  the  close  proximity  of  a  structure  and  they  require  an  extensive  set  of 
constants  which  depend  on  the  charge  material,  which  can  be  very  difficult  to  obtain. 
Another  issue  with  the  similitude  equations  and  incompressible  fluid  analytical  models  is 
the  lack  of  energy  loss  predictions.  Various  analytical  models  have  been  developed  in  the 
past  to  predict  the  time  history  behaviour  of  a  free  field  gas  bubble.  This  report  compares 
the  bubble  radius  time  history  of  a  series  of  analytical  models  to  experimental  fit  results 
reported  by  Swift  and  Decius  [1]. 
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2  Empirical  Models 


2.1  Similitude  Equations 

The  similitude  equation  for  the  first  maximum  bubble  radius  [2]  is  as  shown  in  Eq.  (1). 


R 


max 


'  w  Y3 

v£>  +  loJ 


(1) 


where  Rmax  is  the  maximum  bubble  radius  in  meters,  D  is  the  charge  depth  in  meters,  W  is 
the  charge  weight  in  kilograms,  and  Kb  is  a  constant  depending  on  the  charge  material. 
The  similitude  expression  for  determining  the  first  bubble  period  [2],  Tb,  in  seconds  is 
estimated  by  Eq.  (2). 


Th  =K* 


{D  + 10)% 


(2) 


where  K$  is  a  constant  which  depends  on  the  charge  material.  The  rise  in  the  gas  bubble 
over  its  first  bubble  period,  Ad,  is  given  by  Eq.  (3),  as  show  by  van  Aanhold  [2],  It 
should  be  noted  that  this  is  equation  is  specific  to  TNT. 


Ad  =  12.4 


W 

(D  +  10) 


(3) 


The  volume  of  the  gas  bubble  can  be  reasonably  approximated  by  a  half  sine  function  as 
shown  by  Eq.  (4),  where  V0  is  the  initial  charge  volume  and  Vmax  is  the  maximum  bubble 
volume  calculated  using  the  maximum  bubble  radius  assuming  a  spherical  bubble,  as 
shown  in  Eq.  (4). 


V(t)=V0+{V max-^Jsin 


(  M  ^ 


\Jb) 


4  3 

^"niax  —  "y 
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From  the  volume  time  relation  of  the  bubble  the  radius  time  history  can  be  determined  by 
Eq.  (5)  assuming  a  spherical  gas  bubble. 
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The  vertical  position  of  the  gas  bubble  center,  which  describes  its  vertical  migration,  can 
be  approximated  by  Eq.  (6). 


~(t)  =  -d  +  M 
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2.2  Similitude  Model  With  Energy  Loss 

This  approach  used  the  first  cycle  period  constants  with  an  energy  loss  function  to 
determine  the  radial  and  migration  time  histories.  The  energy  loss  function,  given  by  Eq. 
(7),  was  provided  by  van  Aanhold  [2],  which  was  based  on  reduction  values  used  by 
Geers  and  Hunter  [7]. 

2  =  0.2575  +  (1  -  0.2575)  exp(-  0.8148™ )  (7) 

Here  X  has  an  initial  value  of  1  (/«= 0),  and  m  is  the  cycle  number  starting  at  0  for 
detonation  and  with  m=  1  for  the  first  gas  bubble  minimum.  The  energy  loss  function  is 
implemented  by  reducing  the  charge  mass  at  each  gas  bubble  minimum,  such  that  the 
new  charge  weight  is  XW.  This  new  charge  weight  is  input  into  Eqs.  (1)  through  (3)  to 
determine  the  maximum  gas  bubble  radius,  the  bubble  period,  and  the  migration  during 
the  next  cycle.  This  energy  loss  function  results  in  a  decrease  in  the  maximum  gas 
bubble  radius  and  periods  by  an  order  of  (X  ). 

2.3  Experimental  Fit  Model 

The  average  similitude  constants  for  the  charge  material  used  in  this  study  are  shown  in 
Table  1. 

Table  1:  Similitude  constants  for  TNT  in  MKS  units 


Constant 

TNT 

1st  Cycle 

2nd  Cycle 

3rd  Cycle 

k5 

2.11 

1.57 

1.33 

k6 

3.36 

2.27 

1.76 

The  three  sets  of  constants  for  TNT  are  based  on  an  experimental  study  by  Swift  and 
Decius  [1],  who  experimentally  measured  the  radius  and  period  of  the  gas  bubble.  The 
similitude  constants  were  determined  as  the  values  producing  the  best  fit  to  the 
experimental  data.  For  some  experiments  the  fits  were  produced  for  the  first  three  bubble 
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periods,  but  for  most  experiments  there  were  only  measurements  for  the  first  one  or  two 
bubble  periods.  The  experimental  set-up  and  resulting  similitude  constants  for  all 
experiments  for  TNT  are  shown  in  Annex  A. 
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3  Analytical  Models 


Experimentally  it  is  extremely  difficult  to  monitor  the  position  and  radius  of  the  gas 
bubble  created  by  underwater  explosions.  Therefore,  in  order  to  help  validate  numerical 
approaches  for  solving  the  gas  bubble  behaviour,  analytical  expressions  have  been 
derived  for  predicting  the  time  history  of  the  gas  bubble  radius  and  vertical  position. 
Although  several  approaches  are  available  there  are  limitations  in  the  predictive  models. 
Like  the  similitude  equations,  the  analytical  models  do  not  predict  the  influence  of  a  close 
proximity  target.  The  various  equations  of  motion  (EOM)  compared  in  this  report 
include  the  Rayleigh-Plesset  equation  which  was  original  developed  by  Lamb  [3], 
Herring  EOM  [4],  Kirkwood  and  Bethe  (KB)  EOM  [5],  Keller  and  Kolodner  (KK)  EOM 
[6],  three  EOM  developed  by  Geers  and  Hunter  [7,8,9],  and  an  EOM  developed  at  TNO 
[2]. 

3.1  Lamb  EOM 

The  Lamb  EOM  [3],  shown  in  Eq.  (8),  was  one  of  the  original  analytical  solutions  for 
predicting  the  gas  bubble  radius.  This  approach  assumes  the  bubble  is  stationary  about  its 
center,  such  that  the  gas  bubble  moves  in  the  radial  direction  only  and  has  no  vertical 
migration.  It  assumes  the  fluid  is  incompressible  and  a  free-field  condition. 

PP  +  =  Pw  (pgas  " Pair  +  PwSZ) 


Where  R ,  R  and,  R  are  the  instantaneous  bubble  radius,  radial  velocity,  and  radial 
acceleration  respectively,  pw  is  the  water  density,  pgas  is  the  pressure  in  the  gas  bubble, 
pair  is  the  atmospheric  pressure,  g  is  gravity,  and  z  is  the  depth  of  the  charge. 

Due  to  the  incompressible  fluid  formulation  of  the  Lamb  EOM  there  is  no  energy  loss 
built  into  the  model,  therefore  an  external  energy  loss  function  must  be  applied.  The 
energy  loss  function  is  the  same  factor  as  that  used  for  the  similitude  equations,  with  the 
loss  function  A  being  applied  to  the  internal  gas  pressure  of  the  bubble,  as  shown  by  Eq. 
(7)  in  Section  2.2. 

3.2  Herring  EOM 

The  Herring  EOM  [4],  Eq.  (9),  accounts  for  wave  effects  (i.e.  compressibility)  in  the 
water  but  not  in  the  gas  bubble.  The  compressibility  is  accounted  for  through  the 
addition  and  inclusion  of  a  finite  fluid  sound  speed,  c.  By  setting  the  sound  speed  to 
infinity,  the  Herring  EOM  reduces  to  the  Lamb  EOM  for  an  incompressible  fluid. 
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Pg  is  the  first  time  derivative  of  the  pressure  in  the  gas.  With  the  inclusion  of 
compressibility,  energy  losses  can  occur.  The  most  significant  energy  loss  term  is  the 
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term  containing  P  ;  the  term 
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4R 

3c 


gives  negligible  energy  loss,  but  creates  an 


asymmetry  in  the  radius  time  history,  whereby  the  contraction  of  the  bubble  is  slower 
than  its  growth. 


3.3  Kirkwood  and  Bethe  (KB)  EOM 

Similar  to  the  Herring  EOM,  the  KB  EOM  [5],  Eq.  (10),  accounts  for  wave  effects  in  the 
liquid  and  can  be  implemented  for  either  an  incompressible  or  compressible  fluid.  As 
this  equation  shows  there  is  very  little  difference  between  this  equation  and  that  produced 
by  Herring. 
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The 


1  +  ^ 
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term  on  the  right  hand  side  makes  a  slight  difference  due  to  its  influence  on 


the  damping  term  containing  Pg  .  This  results  from  keeping  the  enthalpy  at  the  bubble 
surface  throughout  the  derivation,  and  not  using  an  approximation  to  the  enthalpy  as 


implemented  in  the  Herring  and  KK  EOM.  The  other  terms  containing 
negligible  effects. 


rk' 


produce 


3.4  Keller  and  Kolodner  (KK)  EOM 

The  KK  EOM  [6],  Eq.  (11),  looks  similar  to  the  KB  model  but  is  very  similar  in 
behaviour  to  that  produced  by  Herring,  as  the  major  damping  term  is  the  same.  The  KK 
EOM  accounts  for  the  fluid  compressibility  by  defining  a  finite  value  for  the  speed  of 
sound  in  the  water. 
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Both  the  Herring  and  KK  EOM  can  be  grouped  into  a  single  one-parameter  equation, 
shown  in  Eq.  (12),  which  was  referred  to  as  the  Keller-Herring  Equation  by  Prosperetti 
and  Lezzi  [10],  where  A  is  an  arbitrary  parameter  with  a  value  between  zero  and  one. 
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Setting  A  to  zero  gives  the  KK  EOM;  setting  A  to  1  gives  the  Herring  EOM. 

3.5  Geers  and  Hunter  EOM  Models 

Over  the  past  decade  and  a  half  Geers  and  Hunters  have  developed  several  models  to 
predict  the  gas  bubble  behaviour  for  underwater  explosions.  For  the  purposes  of  this 
study  three  EOM  variations  were  used.  Two  approaches,  the  DA  and  coupled  DA 
account  for  compressibility  in  the  fluid,  where  the  third  model,  DAA  EOM,  accounts  for 
wave  effects  in  both  the  fluid  and  gas  bubble. 

3.5.1  GH-DAEOM 

For  the  GH-DA  EOM  [7],  only  the  dilatational  motion  of  the  gas  bubble  is  considered,  as 
given  by  Eq.  (13). 
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If  the  DA  EOM  is  used  for  an  incompressible  fluid,  it  is  reduced  to  the  Lamb  EOM. 

3.5.2  Coupled  GH-DA  EOM 

The  coupled  GH-DA  EOM  [8]  is  the  GH-DA  approach  with  migration  motion  coupled  to 
the  bubble  radius  EOM.  The  dilatation  and  migration  EOM  are  shown  in  Eqs.  (14)  and 
(15)  respectively. 
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It  was  found  that  the  inclusion  of  compressibility 

v  ~  "V 

small  or  negligible  effect;  therefore  it  was  not  included  for  this  study. 


in  the  migration  EOM  has  a  very 


3.5.3  GH-DAA  EOM 

The  GH-DAA  EOM  [9],  Eq.  (16),  is  the  only  analytical  expression  that  takes  into  account 
wave  effects  in  both  the  external  fluid  and  the  gas  bubble.  This  model  does  not  account 
for  migration  of  the  gas  bubble. 
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Where  £  =  pgcg / p,c,  is  the  specific-acoustic-impedance  ratio  and  subscripts  /  and  g 

indicate  the  water  and  gas  bubble  respectively.  The  GH-DAA  EOM  reduces  to  that  of 
Keller  and  Kolodner  if  the  wave  effects  in  the  gas  are  not  included.  This  is  achieved  by 
setting  pg  =  pg  =  ^  ^  =  0  in  Eq.  (16).  Including  only  the  wave  effects  in  the  gas  are 

achieved  by  setting  £  =  0,  c,  oo  and  £ ct  =  {pg / pt)cg  in  Eq.  (16).  This  model  cannot  be 

implemented  with  incompressible  fluids  unless  all  wave  effects  are  removed  within  the 
fluid  and  gas  bubble,  which  was  not  considered  in  this  study. 

3.6  Modified  Geers-Hunter  (TNO)  EOM 

Based  on  the  theory  of  the  Geers-Hunter  DA  EOM,  van  Aanhold  [2]  developed  analytical 
EOM  for  the  dilatation  and  migration  of  the  UNDEX  gas  bubble.  These  will  be  referred 
to  as  the  TNO  EOM  within  this  report,  and  are  shown  in  their  general  form  in  Eqs.  (17) 
and  (18)  for  the  dilatation  and  migration  motion  respectively. 
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Compared  to  the  DA  EOM,  additions  by  van  Aanhold  include: 


(18) 
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•  A  kinetic  energy  coupling  switch,  K.  This  couples  the  kinetic  energy  of  the 
dilatation  motion  to  the  migration  motion 

•  Gas  inertia  effect  in  the  dilatation  and  migration  motions,  shown  by  inclusion  of 
the  charge  weight,  W,  in  the  acceleration  terms. 

Similar  to  the  other  models  that  have  been  developed  for  compressible  fluids,  the  TNO 
EOM  can  be  used  for  an  incompressible  fluid  by  setting  c,  -» ® .  The  kinetic  energy 
coupling  between  the  dilatation  and  migration  motions  can  be  turned  off  by  setting  K=  0. 

3.7  Equation  of  State 

The  equation  of  state  is  used  to  calculate  the  pressure  within  the  gas  bubble  throughout 
the  growth  and  collapse  cycles.  It  is  desirable  to  have  an  adiabatic  equation  of  state 
(EOS)  for  the  ease  of  calculations.  Previous  studies  carried  out  by  Jones  and  Miller  (JM 
EOS)  [11]  and  Lee,  Horning,  and  Kury  (JWL  EOS)  [12]  who  extended  the  work  of  Lee 
and  Miller  developed  adiabatic  EOS  for  dense  gas.  Both  of  these  forms  produced 
consistent  models  for  the  adiabatic  expansion  of  the  gas  bubble  for  densities  in  the  range 
of  1.0-1.63  g/cm3  and  at  pressures  below  100  MPa.  For  these  conditions  the  EOS  can  be 
simplified  to  the  ideal-gas  fonn  given  by  Eq.  (19),  which  is  implemented  due  to  its 
simplicity  compared  to  the  more  advanced  JM  and  JWL  EOS. 


II 

VCN 

II 

M 
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3  r 


(19) 


where  P  is  the  pressure  within  the  gas  bubble,  X  is  the  same  energy  loss  function  that  is 
applied  to  the  similitude  plus  energy  loss  models,  K  is  the  adiabatic  charge  constant,  V  is 
the  bubble  volume,  R  is  the  bubble  radius,  y  is  the  constant  ratio  of  specific  heats, 
subscript  g  denotes  the  gas  bubble,  and  subscript  c  denotes  the  initial  charge.  For 
compressible  fluid  models  the  energy  loss  function  is  always  equal  to  1.  For 
incompressible  fluid  models  the  energy  loss  function  is  the  same  as  it  is  for  the  similitude 
models. 

The  adiabatic  charge  constant  and  ratio  of  constant  specific  heats  for  the  incompressible 
fluid  models  can  be  determined  from  either  the  JM  or  JWL  EOS.  Due  to  the  energy  loss 
associated  with  the  compressible  fluid  models,  the  charge  constant  for  compressible  fluid 
models  must  be  significantly  larger  than  that  for  incompressible  fluids.  Similar  to  van 
Aanhold  [2],  the  charge  constant  for  the  compressible  fluid  models  in  this  study  were 
selected  such  that  the  first  bubble  radius  and  period  agree  reasonably  well  with  the  first 
cycle  of  the  similitude  equations.  The  ratio  of  constant  specific  heats  is  the  same  as  that 
used  for  the  incompressible  fluid  models. 
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4  Solution  Methods 


4.1  Similitude  Equations 

The  similitude  equations  provide  three  approaches  for  calculating  the  radial  and 
migration  motion  time  histories  with  Eqs.  (4)  through  (6).  One  uses  these  equations 
without  energy  losses;  the  second  employs  the  energy  loss  function  as  described  in 
Section  2.2. 

The  third  approach  uses  the  experimental  results  published  by  Swift  and  Decius  [1]  to 
directly  calculate  the  radius  and  period  of  each  bubble  cycle,  referred  to  as  the  “Exp.  Fit” 
model.  At  the  first  gas  bubble  minimum  the  similitude  constants  are  changed  to  the 
values  for  the  second  cycle.  The  same  approach  is  used  at  the  second  gas  bubble 
minimum.  For  the  experimental  cases  which  did  not  have  best  fits  for  the  first  three  gas 
bubble  cycles,  the  average  values  shown  in  Table  1  were  used.  This  approach  produces 
results  that  are  very  close  to,  and  can  be  considered  to  be  identical  to,  the  experimental 
measurements. 


4.2  Analytical  Equations  of  Motion 


The  EOM  for  the  analytical  models  in  Eqs.  (8)  through  (18)  are  nonlinear  second-order 
differential  equations  that  need  to  be  solved  by  numerical  integration.  If  the  EOM  can  be 
expressed  in  the  first  order  formj  =  the  numerical  solution  is  obtained  using  a 

fourth-order  Runge-Kutta  method  [13],  which  has  the  basic  form  as  shown  in  Eqs.  (20) 
and  (21). 


T/+i  =y,+A? 
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where  yt  is  the  variable  at  the  start  of  the  time  step,  At  is  the  time  step,  and  j{y,  t)  is  the 
first  order  differential  variable  evaluated  at  y  and  t. 
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**  At  \ 

A+l/2  =  A  +  y  /  v  f+l/2  ’  0+1/2  / 

y*+ 1  =yi+Atf(yM/2’ti+i/2) 


(21) 


A  second  order  differential  equation  (DE)  is  reduced  to  a  series  of  first  order  DEs  by 
creating  a  dummy  variable,  in  this  example  q,  identical  to  the  radial  velocity  of  the  gas 
bubble.  The  dummy  variable  is  substituted  into  the  original  DE.  Eq.  (22)  shows  how  this 
is  done  to  produce  two  first  order  equations  for  the  Lamb  EOM. 
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(22) 
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The  first  order  differential  equations  are  solved  with  the  initial  conditions  to  determine 
Initially  the  gas  bubble  is  at  rest  with  an  initial  radius  equal  to  the  initial  charge 

radius,  such  thatk!  =  Rc,  R=q=q=  0  .  f(ynt,)  is  then  substituted  into  Eq.  (21)  to  solve 

foryI+1/2 .  This  value  is  used  to  update  all  variables  such  as  bubble  pressure.  Eq.  (22)  is 

used  with  the  updated  variables  to  determine  /(y*+i/2>  *1+1/2)-  This  process  is  repeated  using 

i{!^  J  ^:sf:  \ 

Eqs.  (21)  and  (22)  to  determine yi+y2 and f[yi+ i/2>h+i/2/-  Eqs.  (21)  and  (22)  are  used  to 

determine  y*+i  and/(y*+1,h+1).  Eq.  (20)  was  used  to  obtain  yi+i  and  the  process  is 
repeated  for  the  next  time  step.  The  Runge-Kutta  subroutine  for  the  Lamb  EOM  is 
shown  in  Appendix  B.  It  was  found  that  a  time  step  size  of  5.0x10-7  seconds  resulted  in 
convergence,  and  this  was  therefore  used  for  the  numerical  integration. 

For  the  Lamb  EOM  the  same  energy  loss  function  as  that  used  for  the  similitude 
equations  is  applied  during  the  calculations  to  reduce  the  adiabatic  charge  constant  by  the 
factors  given  by  Eq.  (7).  Without  the  use  of  the  energy  loss  function,  there  would  be  no 
reduction  in  the  peak  bubble  radius  throughout  the  time  history.  For  the  compressible 
fluid  models  the  energy  loss  is  inherent  within  the  EOM,  so  there  was  no  additional 
energy  loss  mechanism  considered  for  this  study. 

4.3  Constants 

The  Swift  and  Decius  [1]  report  gives  no  indication  of  the  charge  density  or  any  other 
material  properties  other  than  the  charge  type.  For  this  reason  two  sets  of  charge 
constants  were  used  for  comparison  purposes.  Cole  [14]  uses  a  gas  exponent,  y=  1.25  for 
TNT,  and  therefore  this  was  considered  for  one  property  set.  The  adiabatic  charge 
constant,  K,  was  selected  to  get  reasonable  agreement  between  the  first  bubble  period  and 
maximum  radius  for  the  compressible  and  incompressible  (Lamb  EOM)  fluid  cases.  The 
second  set  of  constants  is  those  used  by  van  Aanhold  [2]  in  his  study.  The  constants  used 
for  the  analytical  models  are  shown  in  Table  2. 
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Table  2:  Constants  used  for  solving  the  analytical  models 


Variables 

Set  1 

Set  2 

Explosive  Density,  pc  (kg/m  ) 

1500 

1630 

Gas  exponent,  y 

1.25 

4/3 

Adiabatic  charge 
constant  (MPa),  K 

Compressible 

1350 

2455 

Incompressible 

740 

1050 

Drag  Coefficient,  Cd 

2.25 

2.25 

Density  of  seawater  (kg/m  ),  pw  or  pi 

1025 

1025 

Sound  velocity  in  water  (m/s),  c 

1500 

1500 

Acceleration  of  gravity  (m/s  ),  g 

9.80665 

9.80665 

Air  pressure  (kPa),  Pair 

101.325 

101.325 

4.4  Code  Development 

The  visualization  software  IDL  was  used  to  develop  a  code  for  calculating  the  gas  bubble 
radius  and  migration,  and  the  time  histories  for  each  of  the  analytical  and  similitude 
approaches.  A  different  subroutine  was  developed  for  each  of  the  similitude  and 
analytical  models.  For  the  purpose  of  this  study  only  the  radial  and  migration 
displacements  have  been  output  from  the  subroutines,  but  could  easily  be  changed  to 
output  all  required  variables. 
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5  Analytical  and  Similitude  Results 


This  section  compares  the  radial  motion  time  history  for  the  different  similitude 
approaches  as  well  as  the  various  analytical  models. 

5.1  Similitude  Energy  Loss  Model  vs.  Experimental  Fit  Model 


Figure  2  shows  a  comparison  of  the  experimental  measurements  and  the  energy  loss 
function  used  in  conjunction  with  the  similitude  equations. 


Figure  2:  Comparison  of  experimental  similitude  response  and  similitude  with  energy  loss 

function  (a)  Case  G5F  (h)  Case  G70F 

From  Figure  2  it  can  be  seen  that  the  energy  loss  function  does  not  provide  the  required 
amount  of  reduction  in  the  bubble  period  or  radius.  The  energy  loss  function  will  result 
in  a  period  and  radius  reduction  of  approximately  16  percent  for  the  second  bubble  cycle 
and  an  additional  10  percent  for  the  third  bubble  cycle.  From  the  experimental  fit  it  was 
found  that  on  average  the  gas  bubble  period  is  reduced  by  approximately  25  percent  and 
another  1 1  percent  over  the  first  and  second  bubble  cycles  respectively.  The  bubble 
radius  over  the  first  and  second  bubble  cycles  are  reduced  by  32  and  15  percent 
respectively. 

5.2  Analytical  Models 

The  following  sections  compare  the  bubble  radius  time  history  for  the  various  analytical 
models,  including  the  effect  of  fluid  compressibility  and  the  charge  properties  on  the 
radius  time  history.  Figure  3(a)  shows  a  typical  gas  bubble  migration  time  history.  Due 
to  the  small  bubble  migration,  less  than  1  percent  of  the  total  depth,  the  coupling  of 
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bubble  depth  and  radius  had  a  negligible  effect  on  the  gas  bubble  radius  time  histories,  as 
shown  by  Figure  3(b). 


(a)  (b) 

Figure  3:  Effect  of  including  bubble  migration  in  the  analytical  model  (a)  Typical  gas  bubble 
migration  (b)  Radius  for  a  migration  +  dilatation  model  and  dilatation  only 


5.3  Effect  of  Fluid  Compressibility 

Figure  4  shows  the  bubble  radius  time  history  of  all  analytical  models  for  the  depth  and 
charge  size  for  test  case  G5F.  Figure  4(a)  shows  that  all  compressible  gas  models  reduce 
to  the  Lamb  EOM  when  the  fluid  becomes  incompressible. 

-  Lamb  -  Lamb 


(a)  (b) 

Figure  4:  Comparison  of  analytical  bubble  radius  time  history  for  a  charge  size  and  detonation 
depth  equivalent  to  test  G17F(a)  Incompressible  Fluid  (b)  Compressible  Fluid 
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Figure  4(b)  shows  a  comparison  of  all  the  analytical  models  with  the  water  considered  as 
a  compressible  fluid,  with  the  exception  of  the  Lamb  EOM  which  does  not  have  any 
compressible  fluid  modelling  capability.  These  results  are  grouped  according  to 
compressibility.  Considering  just  water  compressibility  the  various  models  show  little 
variation  for  a  detonation  depth  of  93  m.  The  only  models  that  show  a  significant 
difference  are  the  Lamb  EOM  which  does  not  consider  water  compressibility  and  the 
Geers-Hunter  DAA  (GH-DAA)  model  which  includes  compressibility  in  the  water  and 
the  gas  bubble.  The  additional  energy  loss  factors  considered  in  the  GH-DAA 
fonnulation  significantly  reduces  the  maximum  gas  bubble  radius  and  period  when 
compared  to  the  other  analytical  approaches. 

5.4  Analytical  Models  vs.  Experimental  Fit  Model 

Figure  5  through  Figure  8  compares  analytical  models  with  experimental  fit  results.  GiF 
indicates  the  use  of  the  similitude  constants  for  the  ith  test  number  from  Swift  and  Decius 
[1].  Exp  Fit  indicates  the  experimental  fit  constants  for  all  three  bubble  cycles  are  used. 
Sim  +  Loss  indicates  the  first  bubble  cycle  similitude  constants  are  used  along  with  the 
energy  loss  function.  For  the  analytical  models  it  is  the  name  of  the  model  followed  by 
the  property  set  from  Table  2,  where  Propl  indicates  property  set  1  and  Prop2  is  for 
property  set  2. 


Figure  5:  Comparison  of  experimental  similitude  response  and  the  Lamb  EOM  with  energy  loss 

function  (a)  93  m  depth  (b)  1 79  m  depth 

From  Figure  5  it  can  be  seen  that  the  general  behaviour  of  the  Lamb  EOM  follows  that  of 
the  similitude  equation,  with  a  slight  improvement  in  the  bubble  period  predictions  for 
the  Lamb  EOM.  The  same  energy  loss  function  is  applied  to  both  of  these  models.  The 
energy  loss  function  does  not  provide  the  required  energy  reduction  for  the  Lamb  EOM 
as  the  peak  radii  and  bubble  periods  of  the  experimental  fits  are  significantly  less. 
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(a)  (b) 

Figure  6:  Comparison  of  experimental  similitude  response  and  the  GH-DA  EOM  (a)  93  m  depth 

(b)  1 79  m  depth 


(a) 


(b) 


Figure  7:  Comparison  of  experimental  similitude  response  and  the  TNO  EOM  (a)  93  m  depth  (b) 

179  m  depth 


Comparing  Figure  6  and  Figure  7  shows  that  including  the  gas  inertia  effects  in  the  TNO 
EOM  results  in  negligible  differences  in  the  gas  bubble  behaviour  at  the  depths 
considered.  These  figures  also  show  that  the  energy  loss  inherent  within  these  models  is 
not  sufficient  compared  to  the  similitude  plus  energy  loss  or  experimentally  fit  results. 
The  compressible  fluid  models  over  predict  the  radii  by  27  and  54  percent  for  the  second 
and  third  gas  bubble  maximums  and  the  periods  by  11  and  19  percent  compared  to  the 
experimental  fit  results. 
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Figure  8:  Comparison  of  experimental  similitude  response  and  the  GH-DAA  EOM  (a)  93  m  depth 

(b)  1 79  m  depth 

Figure  8  shows  that  the  energy  loss  inherent  within  the  DAA  EOM  formulation  produces 
slightly  better  bubble  radius  results  than  that  shown  by  the  Lamb  EOM  and  significantly 
better  results  than  the  other  compressible  fluid  models.  The  bubble  period  predicted  by 
the  DAA  EOM  is  not  as  large  as  that  imposed  by  the  energy  loss  function  however  it  is 
larger  than  that  produced  by  the  other  compressible  fluid  analytical  models. 

In  general  the  two  property  sets  have  only  a  minor  effect  on  the  later  bubble  cycles,  as 
shown  by  the  convergence  of  the  time  histories  during  the  third  bubble  cycle. 
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6  Conclusion 


An  extensive  study  on  the  analytical  equations  of  motion  available  for  predicting  the  gas 
bubble  behaviour  was  conducted  to  determine  which  models  and  assumptions  provide 
better  predictions.  Using  the  visualization  software  IDL  the  various  analytical  models 
and  similitude  equations,  a  code  was  developed.  This  code  outputs  the  radial  and 
migration  time  history  for  the  various  equations  of  motion  and  similitude  approaches. 

This  study  has  shown  that  although  the  general  gas  bubble  behaviour  is  captured  by  each 
of  the  analytical  models  there  are  still  significant  gaps  in  the  predictions.  It  was  found 
that  the  energy  loss  function  proposed  by  van  Aanhold  [2]  does  not  adequately  predict  the 
actual  loss  that  occurs  in  the  experiments  when  used  in  conjunction  with  the 
incompressible  fluid  models.  It  is  also  apparent  that  the  energy  loss  inherent  in  the 
compressible  fluid  models  does  not  fully  predict  the  energy  loss  observed  in  experiments. 
The  percent  error  in  the  gas  bubble  period  and  maximum  bubble  radii,  compared  to  the 
experimental  fits  for  a  detonation  depth  of  approximately  93  m  and  over  150  m,  are 
shown  in  Table  3  and  Table  4  respectively.  The  compressible  EOM  columns  are  the 
average  of  the  various  models  that  consider  only  compressibility  in  the  surrounding  fluid. 

Table  3:  Percent  error  in  the  gas  bubble  periods  and  maximum  radii  predicted  by  the  analytical 
models  compared  to  the  experimental  fits  for  a  detonation  depth  of  approximately  93  m 


Cycle 

Number 

Lamb  EOM  (%) 

Compressible  EOM  (%) 

GH-DAA 

EOM  (%) 

Period 

Radius 

Period 

Radius 

Period 

Radius 

1st 

1.32 

1.64 

1.01 

1.60 

2.39 

6.59 

2nd 

3.10 

22.3 

10.9 

27.1 

6.78 

15.4 

3rd 

6.47 

41.5 

19.6 

54.5 

15.2 

37.2 

Table  4:  Percent  error  in  the  gas  bubble  periods  and  maximum  radii  predicted  by  the  analytical 
models  compared  to  the  experimental  fits  for  a  detonation  depth  of  over  150  m 


Cycle 

Number 

Lamb  EOM  (%) 

Compressible  EOM  (%) 

GH-DAA 

EOM  (%) 

Period 

Radius 

Period 

Radius 

Period 

Radius 

1st 

1.01 

1.19 

1.86 

1.75 

1.24 

7.01 

2nd 

5.33 

11.3 

14.2 

14.3 

10.6 

3.55 

3rd 

8.13 

39.7 

22.8 

50.2 

19.1 

33.2 

As  these  numbers  show  the  Lamb  EOM  predicts  the  bubble  period  closer  to  the 
experimental  fits  than  the  other  models.  This  shows  that  the  energy  loss  function 
provides  a  larger  loss  than  those  inherent  within  the  compressible  fluid  and  GH-DAA 
EOM.  As  would  be  expected  the  GH-DAA  model  which  accounts  for  fluid  and  gas 
bubble  compressibility  shows  a  significantly  larger  energy  loss  than  the  EOM  considering 
compressibility  in  the  fluid. 

Table  3  and  Table  4  show  that  there  is  a  requirement  to  include  an  energy  loss  function  in 
conjunction  with  the  compressible  fluid  models.  All  models  that  include  only  the 
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compressibility  in  the  fluid  produce  results  very  similar  to  one  another;  therefore  a  single 
energy  loss  function  could  be  developed  for  all  of  these  models.  However,  further  studies 
would  have  to  be  conducted  to  determine  if  the  energy  loss  is  consistent  for  a  range  of 
detonation  depths  and  charge  types  and  sizes.  The  detonation  depth  is  expected  to  be  a 
significant  factor  for  models  that  include  coupling  between  radial  motion  and  migration 
for  shallow  charges. 

A  concern  of  the  approach  used  in  this  report  is  the  method  of  determining  the  adiabatic 
charge  constant  for  the  compressible  fluid  models.  The  charge  constant  required  for  a 
compressible  fluid  is  much  larger  than  that  used  for  incompressible  models  in  order  to  get 
reasonable  agreement  for  the  first  gas  bubble  period  and  maximum  radius.  This  is  due  to 
the  energy  losses  inherent  in  the  compressible  fluid  models.  The  concern  is  that  the  gas 
constant  is  taken  as  a  value  that  fits  the  first  cycle  of  the  gas  bubble.  A  method  that  is 
proposed  by  Geers  and  Hunter  [9]  involves  starting  from  a  time  shortly  after  detonation, 
and  involves  determining  the  bubble  growth  during  the  shock-wave  phase.  This  approach 
allows  for  the  use  of  the  gas  constant  determined  from  the  equations  of  state  to  be  used 
for  the  compressible  gas  models.  This  approach  will  be  introduced  into  the  models  in  the 
future. 
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Annex  A  Swift  and  Decius  Experimental  Results 


The  following  tables  show  the  gas  bubble  periods  and  maximum  gas  bubble  radius  for  the 
experiments  conducted  by  Swift  and  Decius.  Where  the  values  in  the  shaded  cells  were 
calculated  based  on  the  average  Kf  value,  where  i  indicates  the  cycle  number. 


Table  5:  Gas  bubble  period  from  experiments  conducted  by  Swift  and  Decius  [1J. 


Test 

Case 

Charge 

Weight 

(lb) 

Depth 

(ft) 

First 

Period 

(msec) 

ki5 

Second 

Period 

(msec) 

K25 

Third 

Period 

(msec) 

K35 

GIF 

0.651 

343 

27.4 

4.43 

20.1 

3.25 

17.0 

2.75 

G2F 

0.660 

304 

29.8 

4.38 

22.1 

3.25 

18.7 

2.75 

G4F 

0.660 

304 

30.0 

4.41 

22.1 

3.25 

18.7 

2.75 

G5F 

0.662 

305 

29.9 

4.40 

22.1 

3.25 

18.7 

2.75 

G6F 

0.658 

304 

29.9 

4.41 

22.1 

3.25 

18.7 

2.75 

G7F 

0.669 

298 

30.2 

4.36 

22.6 

3.26 

19.1 

2.75 

G8F 

0.663 

304 

29.8 

4.38 

22.7 

3.33 

18.8 

2.75 

G9F 

0.651 

304 

29.6 

4.38 

22.3 

3.30 

18.7 

2.75 

G17F 

0.660 

305 

29.6 

4.38 

23.0 

3.40 

18.7 

2.75 

G18F 

0.660 

302 

29.6 

4.32 

21.6 

3.16 

19 

2.77 

G20F 

0.660 

538 

19.6 

4.47 

13.9 

3.16 

12.1 

2.75 

G21F 

0.651 

593 

19.1 

4.73 

14.1 

3.49 

11.9 

2.94 

G23F 

0.658 

567 

18.1 

4.32 

13.3 

3.17 

10.7 

2.54 

G70F 

0.660 

503 

19.9 

4.30 

15.0 

3.25 

12.7 

2.75 

G71F 

0.658 

463 

21.0 

4.26 

14.9 

3.01 

13.6 

2.75 

G72F 

0.660 

586 

17.9 

4.36 

13.0 

3.17 

11.3 

2.75 

G73F 

0.655 

576 

18.2 

4.38 

13.5 

3.25 

11.4 

2.75 

G74F 

0.660 

556 

18.1 

4.24 

13.9 

3.25 

11.8 

2.75 

G76F 

0.660 

587 

17.7 

4.33 

13.3 

3.25 

11.3 

2.75 

Average  Values 

4.36 

3.25 

2.75 
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Table  6:  Maximum  gas  bubble  radius  from  experiments  conducted  by  Swift  and  Decius[l]. 


Test 

Case 

Charge 

Weight 

(lb) 

Depth 

(ft) 

First 

Rmax  (in) 

ki6 

Second 

Rmax 

(in) 

K26 

Third 

Rmax 

(in) 

K36 

GIF 

0.651 

343 

18.1 

12.6 

12.2 

8.50 

9.50 

6.60 

G2F 

0.660 

304 

18.8 

12.5 

12.8 

8.50 

9.90 

6.60 

G4F 

0.660 

304 

18.8 

12.5 

12.8 

8.50 

9.90 

6.60 

G5F 

0.662 

305 

18.8 

12.5 

12.8 

8.50 

9.90 

6.60 

G6F 

0.658 

304 

18.8 

12.5 

13.2 

8.80 

9.90 

6.60 

G7F 

0.669 

298 

18.9 

12.4 

12.7 

8.40 

10.0 

6.60 

G8F 

0.663 

304 

18.9 

12.5 

12.8 

8.50 

9.90 

6.60 

G9F 

0.651 

304 

18.7 

12.5 

11.4 

7.60 

10.0 

6.70 

G17F 

0.660 

305 

19.4 

12.9 

13.4 

8.90 

9.90 

6.60 

G18F 

0.660 

302 

19.4 

12.8 

12.1 

8.00 

10.6 

7.00 

G20F 

0.660 

538 

15.9 

12.6 

10.6 

8.40 

8.20 

6.50 

G21F 

0.651 

593 

15.7 

12.9 

9.90 

8.10 

7.90 

6.50 

G23F 

0.658 

567 

15.6 

12.6 

9.90 

8.00 

7.80 

6.30 

G70F 

0.660 

503 

16.2 

12.6 

10.9 

8.50 

8.50 

6.60 

G71F 

0.658 

463 

16.6 

12.6 

11.0 

8.30 

8.70 

6.60 

G72F 

0.660 

586 

15.4 

12.5 

11.6 

9.50 

8.10 

6.60 

G73F 

0.655 

576 

15.6 

12.7 

10.5 

8.50 

8.10 

6.60 

G74F 

0.660 

556 

15.7 

12.6 

10.6 

8.50 

8.20 

6.60 

G76F 

0.660 

587 

15.1 

12.3 

11.6 

9.50 

8.10 

6.60 

Average  Values 

12.6 

8.50 

6.60 
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Annex  B  Runge-Kutta  Code  for  the  Lamb  EOM 


Step  1)  Reduced  differential  equation  to  series  of  first  order  differential  equations 


LAMB  EOM :  RR  + 1 R 2  =  pj  (pgas  -  pair  +  pwgz) 
reduces  to : 


f\  =  q  =  — 
R 


-G 


pw  \p 


gas 


\  3  2 
Pair+PWSM--q 


f2  —  R  —  q 


Step  2)  Calculate  f{yl ,  t( )  using  initial  conditions  t  =  0,  q  =  q  =  0  ,  R  =  Rc : 


fl  1  =  calculate _fl_LAMB_EOM(a,  rho_w,  Pg,  gamma,  Rc,  P_amb,  q,  K,  scale) 
f2 1  =  calculate  J2_LA  MB_EOM(q ) 


where  fll  and  f21  are  fl  and  f2  at  the  initial  conditions  respectively.  The 
calculate Jl_LAMB_EOMQ  is  the  call  to  the  subroutines  that  calculate  the  functions// 
and  f2.  The  variables  in  the  brackets  are  values  passed  into  the  subroutines. 

Step  3)  Calculate  the  yM/2  values  using  /(_v,  ,t,  ): 

ql  1  =  calculate _var( q,  dt/2.0,  fl  1) 
al  1  =  calculate _var{ a,  dt/2.0,  f21) 

$ 

where  ql  1  and  all  are  the  yi+ y2  for  fl  and  /2  respectively,  qll  and  all  are  the  updated 

radial  velocity  and  radius  of  the  gas  bubble  respectively,  calculate _var()  is  the 

subroutine  developed  to  calculate  the  variables. 

Step  4)  Calculate  /(y*+i/2 >*1+1/2)  using  Step  3  results 

fl2  =  calculate _fl _LAMB_EOM(all,  rho_w,  Pg,  gamma,  Ro,  P_amb,  ql  1,  K,  scale) 
m  =  calculate J2_LAMB_EOM(q  I  1 ) 

afcjfc  I  *  \ 

Step  5)  Calculate  the  yi+y2  values  using  f\yi+\/2di+\/2  )'■ 

ql2  =  calculate _var{ q,  dt/2.0,  fl2) 
al2  =  calculate _var(a,  dt/2.0,  f22) 

Step  6)  Calculate  f[yZy2’{i+ 1/2)  using  Step  5  values  for  q  and  a 

fl3  =  calculate _fl _LAMB _EOM{a\2,  rho_w,  Pg,  gamma,  Ro,  P_amb,  ql2,  K,  scale) 
f23  =  calculate  J2_LAMB_EOM(ql2 ) 

*  /  **  \ 

Step  7)  Calculate  the  yM  values  using  f\yi+i/2Ai+i/2  )'■ 
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ql3  =  calculate j^ar{ q,  dt,  fl 3) 
al3  =  calculate_var{ a,  dt,  f23) 

Step  8)  Calculate  f[y*+l,ti+])  using  Step  7  values  for  q  and  a: 

fl4  =  calculate J1LA MB_EOM( a  1 3 ,  rho_w,  Pg,  gamma,  Ro,  P_amb,  ql3,  K,  scale) 
£24  =  calculate  J2_LAMB_EOM(ql  13) 

Step  9)  Calculate  the  variables  at  the  end  of  the  time  step  using  Eq.  as: 


9m  =  +  At 

oM  =yi+At 


i/ll  +  i/12  +  i/13  +  i/14 

o  3  3  o 

— /2 1  +  -/22  +  —  /23  +  —  /24 
6  3  3  6 


-  A+i 

=  ^/+i 


Step  10)  Repeat  steps  2-9  for  each  time  step. 
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List  of  symbols/abbreviations/acronyms/initialisms 


ac  orRc 
c 

Cd 
CFD 
D  or  d 
DE 
DND 
DR  DC 

DRDKIM 

EOM 

EOS 

ft 

g 

JM 

JWL 

K 

K5 

k6 

KB 

Kg 

KK 

lb 

m 

MFP 

Pair 

Pg 

4 

R  or  Rg 


Charge  radius 

Speed  of  sound  in  water 

Drag  coefficient 

Computational  fluid  dynamics 

Depth  of  charge  detonation 

Differential  equation 

Department  of  National  Defence 

Defence  Research  &  Development  Canada 
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